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Abstract 

The transport and chemical reactions of solutes are modelled as a cellular automa- 
ton in which molecules of different species perform a random walk on a regular lattice 
and react according to a local probabilistic rule. The model describes advection and 
diffusion in a simple way, and as no restriction is placed on the number of particles at 
a lattice site, it is also able to describe a wide variety of chemical reactions. Assuming 
molecular chaos and a smooth density function, we obtain the standard reaction- 
transport equations in the continuum limit. Simulations on one- and two-dimensional 
lattices show that the discrete model can be used to approximate the solutions of 
the continuum equations. We discuss discrepancies which arise from correlations be- 
tween molecules and how these disappear as the continuum limit is approached. Of 
particular interest are simulations displaying long-time behaviour which depends on 
long- wavelength statistical fluctuations not accounted for by the standard equations. 
The model is applied to the reactions a + b ^ c and a + b c with homogeneous 
and inhomogeneous initial conditions as well as to systems subject to autocatalytic 
reactions and displaying spontaneous formation of spatial concentration patterns. 
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1 Introduction 



The transport of aqueous solutions of contaminants in geological media is inextricably 
coupled with a rich variety of physical and chemical processes. The solutes may decay 
radioactively, react chemically with each other, sorb on solid surfaces or change the poros- 
ity of the host rock by precipitation/dissolution. Modelling an aqueous system in its full 
complexity constitutes, conceptually and mathematically, a formidable task. A model is 
essentially defined by selecting, on grounds of usefulness and economy, an appropriate set 
of dependent variables (e.g. aqueous species concentrations) and writing down the laws 
(e.g. mass conservation, law of mass action) these variables satisfy. A model intended 
for practical applications has to be translated subsequently into a usable and efficient 
computer code. 

The model presented here takes advantage of the fact that the migration and chemical 
transformation of aqueous species consists, at the microscopic level, of processes taking 
place in parallel (i.e. simultaneously at many locations) and locally (i.e. involving molecules 
within a small spatial neighbourhood). Each species is represented by a (large) number 
of 'particles' that move and react according to simple rules mimicking the microscopic 
behaviour of the actual molecules. The model is implemented as simple computer code 
suitable for massively parallel computers. Our approach contrasts with the standard mod- 
elling approach which solves systems of non- linear partial differential equations (PDE's) 
for macroscopic variables. The latter approach neglects microscopic details by effectively 
averaging solute properties over a small macroscopic volume and dealing only with their 
local mean values. 

Reaction-transport processes are modelled here, following Refs. ^, as a cellular 
automaton ( CA ). In general, a CA is a dynamical system which consists of a discrete- 
valued field defined at the sites of a regular lattice and evolving in discrete time steps, 
with the field value at a site being determined for the next time step by its present values 
in a neighbourhood of the site of interest [^]. In this work, particles reside on the sites 
of a regular lattice; they move randomly between neighbouring sites, in discrete time 
steps, and react with a certain probability upon meeting. Formally, this is described 
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by a set of occupation numbers, i.e. an integer-valued field with a species label, giving 
the number of particles of the different species at each lattice site and evolving in time 
according to a local rule. The evolution rule consists of a 'transport' operation followed 
by a 'reaction' operation. The two operations are repeated iteratively for consecutive time 
steps. Each operation amounts to applying a simple algorithm on all sites of the lattice 
simultaneously. The 'transport' algorithm prescribes the probability with which a particle 
will move to a neighbouring site; the probability may depend on species label, location on 
the lattice, spatial direction or time; in the end, each particle performs a random walk. 
The 'reaction' algorithm, on the other hand, provides the probability with which a given 
combination of occupation numbers will lead to chemical reaction. We define macroscopic 
quantities by averaging over an ensemble of independent copies of the system of interest. 
Particle density, defined as an ensemble average of the occupation number, obeys a discrete 
evolution equation. If chemical reactions are present in this equation, occupation numbers 
can be entirely eliminated in favour of locally averaged particle densities only if certain 
important conditions are met. As will be discussed in detail in Section ^, one has to 
assume that the chemical reactions do not give rise to correlations among different species 
and that the spatial dependence of macroscopic quantities is sufficiently smooth. The 
continuum limit of the resulting equation (i.e. when lattice spacing and time step go to 
zero) contains the standard advection, diffusion and reaction terms of the macroscopic 
PDE's 0. The probabilities of motion and reaction can be chosen so that the desired 
macroscopic physical and chemical parameters are obtained. 

The standard way to model reaction-transport phenomena consists in applying (i) local 
mass conservation to derive a set of partial differential equations, and (ii) local chemical 
equilibrium, according to the law of mass action, to derive a set of non-linear algebraic 
equations (assuming implicitly that chemical equilibrium is attained instantaneously across 

the averaging volume on the time scale of the transport processes). The ensuing system 

^For brevity we shall talk about 'difltusion', although, in the cases of interest to us, it is mechanical dis- 
persion, rather than molecular diffusion, that accounts overwhelmingly for the coefficient of hydrodynamic 
dispersion. To the extent that mechanical dispersion can be described by an effective diffusion term, the 
two processes are indistinguishable from the modelling point of view; naturally, the interpretation of the 
effective coefficient is very different in the two cases p|. We note incidentally that we find it sufficient for 
our purposes to use the term 'species' in a generic sense, although one distinguishes in principle between 
elementary components and composite species. 
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of coupled, non-linear PDE's can be solved by a variety of numerical techniques. A 
common problem is the difficulty in establishing general criteria guaranteeing the stability 
of the numerical algorithms for solving systems of non-linear PDE's; these questions are 
handled on a case by case basis. Concerning the computer implementations of these solvers, 
the memory and time requirements present an extremely difficult task for practical two- 
and three-dimensional problems with present-day computers. As suggested in Ref. 
algorithms that apply transport and chemical reactions sequentially appear to hold more 
promise. We have already seen that this iterative aspect is shared by our model. 

The approach proposed in this work has been motivated by the following arguments: 

(a) We model physicochemical processes at a level intermediate between the macroscopic 
level described by PDE's and the microscopic level of molecular dynamics. The 
'particles' of our model are mathematical abstractions of the actual molecules of 
various chemical species. The number of particles is large enough to make statistical 
concepts meaningful, but is still many orders of magnitude smaller than the real 
number of molecules involved. The CA is intended to describe macroscopic behaviour 
which does not depend on the details of the microscopic dynamics, but follows from 
general properties of the latter. The CA guarantees these essential properties by 
applying them directly to its elementary constituents, in a much more intuitive way 
than the continuum approach. For example, the constraints of mass and momentum 
conservation, together with some smoothness assumptions about the macroscopic 
variables, lead to the Navier-Stokes equations of fluid dynamics i. Similarly, the 
diffusion equation follows essentially from the random nature of collisions at the 
microscopic level. The resulting computer code is simple and allows easy variation 
of the dynamics and the boundary conditions. Moreover, statistical fluctuations are 
inherently present in the CA, whereas PDE's rely on being able to define a physical 
quantity, e.g. particle density, as a continuous variable and thus neglect any local 
variations due to the microscopic nature of the process. This difference is crucial 

when the fluctuations have significant macroscopic consequences, as for example in 

^Smoothness assumptions are indispensable for the derivation of macroscopic PDE's, such as the Navier- 
Stokes equations, from microscopic equations, such as the Boltzmann equation. 
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the reaction a + 6 — > nothing\^^. 

(b) By working with integers, the CA simulation is free of round- off errors, which arise 
from the representation of real numbers by finite computer words and lead to possible 
instabilities of the numerical algorithms used to solve PDE's. This is a significant 
advantage when non-linear reaction terms are present, making the stability of PDE 
solvers hard to establish. All quantities in a CA simulation are intrinsically finite 
and infinities can only arise by extrapolation to some limit (e.g. infinite volume), 
in very much the way this is done with real systems The continuity required 
of the solutions of PDE's is, by contrast, physically inaccurate, as it presupposes 
an infinitesimal limiting procedure (e.g. averaging over a volume that shrinks to a 
point) which does not correspond to reality below a minimal scale (e.g. the volume 
per particle when evaluating particle density). 

(c) CA algorithms are naturally parallel, i.e. they process a large number of data si- 
multaneously. This makes them suitable for massively parallel computers. In par- 
ticular, the simple nature of the elementary physicochemical processes should allow 
simulation on massively parallel computers with relatively unsophisticated proces- 
sors. Moreover, the local character of these processes implies that only a minimal 
amount of communication between physically neighbouring processors will be neces- 
sary. Computers with large numbers of processors (up to tens of thousands) operat- 
ing in parallel and favouring next-neighbour communication are among those leading 
the push towards terafiop performance (10^^ fioating point operations per second) 0]. 
Supercomputers are also being used to solve PDE's with standard numerical meth- 
ods, and the question arises if CA can provide a more efficient alternative (when 
fiuctuation effects are uninteresting). In this respect, we note that the stochastic 
nature of the CA approach makes it in general slower than deterministic methods 
of solving PDE's, if the same time and space discretisation is used. On the other 
hand, these deterministic methods suffer invariably from stability problems, which 
can be overcome, but usually at great expense to the computation time. Here the 
CA approach has the important advantage that it is inherently stable. We there- 
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fore expect that the answer to the above question will depend highly on the specific 
application of interest. 

The systems simulated in this work are conventionally described by reaction-transport 
equations of the form derived in Section ^ (Eq. ( |2.45| )). As discussed above, these equations 
are based on several oversimplifying assumptions about the behaviour of real systems, 
but we shall use them here anyway as a reference point in testing our simulations. We 
justify this choice by the fact that reaction-transport PDE's lie invariably at the basis 
of all known models for the systems in which we are interested. We remain naturally 
open to the possibility that discrepancies between the CA approach and the solutions of 
the reaction-transport equations may arise from some element of microscopic reality (e.g. 
microscopic fluctuations) that the continuum approach fails to capture. This may lead to 
interesting corrections to the standard point of view. When no such corrections arise, the 
CA may be used to approximate the solutions of the continuum equations; it is of course 
essential to study how the result of the simulation converges to the continuous function it 
is approximating, as we refine the discretisation of space and time. The merits of the CA 
approach will then be judged in comparison with other numerical methods. The aim of this 
work is to establish that the proposed CA can simulate a wide variety of physicochemical 
phenomena, ranging from simple annihilation reactions to complex autocatalytic reaction 
schemes leading to pattern formation. The discrete approach will be shown to be capable 
of approximating the solution of the reaction-transport PDE's. In some cases the results 
of the discrete and continuum approaches will disagree and in one of these the discrepancy 
will point to a fundamental shortcoming of the continuum approach. 

CA simulations are currently being performed by scientists in various disciplines, who 
seek in the simplicity of their elementary algorithms the unifying principles behind the 
often complex phenomena they observe Thus, fluid motion is being modelled via 
lattice-gas automata (LGA). Lattice-gas models have existed for several years and they 
were known to display certain hydrodynamic properties but the recent increase in 
their popularity ||l^ has been largely brought about (i) by the advent of powerful super- 
computers, and (ii) by the realisation that they can approximate hydrodynamic PDE's 
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LGA are special cases of deterministic CA. Our approach is more closely related with 
probabilistic LGA models of reaction-diffusion processes 13 1. All LGA models impose 
a limitation on the number of particles per lattice site, in the form of an exclusion principle. 
This is necessary when simulating CA on special purpose computers, which can support a 
small number of bits per site (e.g. CAM-6, [p!4|), and is useful for optimisation purposes on 
vector supercomputers. With computers becoming available that perform floating-point 
operations on many processors in parallel, we see no reason to maintain such a restric- 
tion, in particular since the enforced small number of particles per site severely limits the 
statistics of the simulation and makes either a bigger lattice or more simulations neces- 
sary. In addition, probabilistic models of diffusion with an exclusion principle are not able 



to model advection without introducing at the same time unwanted non-linearities |15]. 
In our model, by contrast, advection emerges naturally from the microscopic rule. An 
exclusion principle also makes it impossible to model local chemical reactions of arbitrary 
complexity. Increasingly complex reactions can be brought into our model with minor 
programming effort. Any model aiming at the description of the transport and chemistry 
of complex solute systems in geological media must be able to incorporate advection and 
arbitrary chemical reactions naturally. The model proposed here appears, in this light, 
most suitable for applications to systems of this kind. 

The paper is organised as follows: In Section ^ the discrete model is introduced and the 
continuum limit is derived with special emphasis on the physical assumptions that lead 
to the desired continuum equations. Alternative microscopic rules and boundary condi- 
tions are also discussed. Extensive simulations of physicochemical systems and detailed 
comparisons with the corresponding differential equations are undertaken in Section ^. 
Discrepancies with the continuum approach will receive particular attention. Finally, in 
Section ^, we discuss computational feasibility, summarise our conclusions and set tasks 
for the future. 
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2 The Model 
2.1 Time Evolution 

Our aim in this section is to formulate a microscopic model for the time evolution of a 
system of particles moving randomly (with a possible bias in a given direction) and reacting 
chemically among themselves. The particles can be thought of as molecules in solution, 
migrating in a porous medium. The model will be general enough to simultaneously 
describe the diffusion and advection of several chemical species with different transport 
properties and subject to various chemical reactions. 

In our model, space and time are discrete: particles reside on the sites of a regular 
lattice, with spacing A, and the system evolves by a sequence of instantaneous transitions 
separated by time r. In each transition, the system moves to a new state according to 
a local rule applied to all sites of the lattice, and the procedure is iterated in discrete 
time steps. Below, we shall define the local rule and proceed to show that, under certain 
conditions, the equations describing the time evolution of the discrete system go over to 
a set of differential equations in the limit X,t ^ 0. 

We define an evolution step to consist of a transport step and a chemical reaction 
step. During the transport step, we consider each lattice site independently and, for each 
particle present, we make a random decision whether it will move or stay stationary, 
with probabilities chosen so as to obtain the desired macroscopic parameters. Having 
redistributed particles this way, we proceed to the chemical reaction step, in which we 
decide for each site independently whether the particles occupying it will react, with 
probabilities reflecting the macroscopic reaction rates. The reaction probability can be 
chosen according to a variety of rules, some of which may give rise to the same macroscopic 
behaviour. A minimal prerequisite for a reaction to take place is that there are sufficient 
particles of each reactant. For the reactions that do go ahead, we remove particles of 
the reactant species and add particles of the product species, according to the reaction 
equation. Having completed the chemical reaction step we proceed to the next evolution 
step and so on. 

We begin the precise formulation of our model by considering a system of particles 
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belonging to various species Sq, (for most purposes the index a will be used interchangeably 
with the name Sa of the species). The particles are located on the sites of a lattice C. 
For simplicity, we consider a one-dimensional lattice. Extensions to higher dimensions are 
straightforward. We focus on a typical site at position x (we shall refer to it for brevity as 
'site x') in the interior of the lattice. At time t, there are Nci{x,t) particles of species a at 
X. We refer to Na{x, t) as the occupation number. The values of Na{x, t) for all species a 
and lattice sites x define the state of the system of particles at time t. We assume a given 
distribution, Na{x,0), of particles at t = 0, which provides the initial condition to our 
problem, and defer to a later point the discussion of boundary conditions. The evolution 
during a time step r, from time t to t + t, amounts to applying on Na{x,t) the product 
C o T of the operators describing transport, T, and chemical reactions, C: 

Nc,{x,t + T) = CoTNa{x,t). (2.1) 

The purpose of the transport operation is to move particles by one lattice spacing to 
the right with probability p, or to the left with probability q {p+q < 1). This is completed 
in two consecutive stages. Each stage consists of a local operation applied simultaneously 
to all lattice sites. Since the action of T is defined for any set of occupation numbers, 
independent of a particular state of the system, we introduce a generic set of occupation 
numbers {J\fa{x) : x G C — dC}, where dC represents the lattice boundary. First, for every 
lattice site x and every species a, a random triplet ^^.^^ = {ix,n\£!x}i-,^x\n^) is drawn for 
each of the Ma{x) particles of species a present at x (n = 1, ...,J\fa{x)). takes one of 
the values (1,0,0), (0, 1,0) or (0,0,1) with probabilities p, 1 — p — g and g respectively. All 
^'s are drawn anew at each update, but an explicit time index is omitted for simplicity. 
These triplets are stored for the second stage. We begin the second stage with a new, 
empty lattice and build the new occupation numbers according to 

1 oo 

TXaix) = E E ^%x,J(.^a{x + jX) - n), (2.2) 

j=-ln=l 

where j, n are integers and ^x+jX n' ^x+jX n' ^x+jX n components of the random 

triplet ^x+jX,n- '^^^ ^-function is defined by 6{m) = 1 for m > 0, and 6{m) = otherwise. 

To show how the empty lattice is filled, we recall that TJ\fa{x) = at the beginning 
of the second stage. The site x is fed with particles from its immediate neighbourhood, 
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i.e. the sites preceding and following x, as well as the site itself (in 2 dimensions, the 
immediate neighbourhood of a site consists of 5 or 7 sites, for a square or triangular 
lattice respectively). On the RHS of Eq. ( |2.2| ), j runs over the sites of this neighbourhood, 
which are x — X,x and x + A. For each one of these sites and for each particle present at 
the same site of the original lattice, we look up the corresponding random triplet stored 
previously. If the component Cx+jXn ™^ particle is added to TJ\fa{x), if it is 0, 
TJ\fa{x) is left unmodified. In physical terms this corresponds to particles moving to the 
right or left or remaining stationary with probabilities p, q and \ —p — q respectively. The 
fact that one component of „ is always 1 and the rest implies that the transport 
operation will place exactly one particle in the immediate neighbourhood of a site for each 
of the particles originally at that site, thus conserving particle number. Then Eq. ( |2.2[ ), 
with Na{x, t) substituted for Ma{x), can be interpreted as yielding the number of particles 
of species a at x, after the transport phase of the update following time t: particles from 
the left (i.e. on site x — A at time t) move to x with probability p, particles from the right 
(i.e. on site x + A at time t) move to x with probability q and particles at x remain there 
with probability 1 — p — q. The last probability can be related to a retention factor in the 
movement of the species a (see Section p. 

Next we define the chemical reaction operator C using the equation: 

R 

CMa{x) = Ma{x) + (^ii^ - ^it^) (2.3) 
r=l 

Here the summation runs over all chemical reactions in the problem at hand and r]x,r is 
a Boolean random variable. Val and Var are the stoichiometric coefficients referring to 
initial and final species of the reaction r respectively. Naturally one can talk about 'initial' 
and 'final' species only if r is a one-way reaction. In fact, most of the reactions we expect 
to encounter are reversible, of the type 

^VaSa ^ YlJ-aSoi, (2.4) 
a K2 

where Ki,K2 are the rate constants and the stoichiometric coefficients are non- 

negative integers. The summations run over all species, but the stoichiometric coefficients 
vanish for species not involved in the reaction. It is convenient to express the above 
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reaction as two one-way reactions, namely 

^gSg > fJ'aSg, (2.5a) 

a a 

^fJ-aSa ^'^aSa- (2.5b) 

a a 

In Eq. (|2.3D R is the total number of one-way reactions thus obtained and the index r 
runs over all such reactions (if there are irreversible reactions, they are added to the list 
as they are). Thus, if the first of the above two reactions occupies the r-th position on 
the list of one-way reactions, we identify Ua and /u^ with Uar and Uar respectively. In 



Eq. ( |2.3| ) rjx^r determines whether reaction r takes place {r]x,r = 1) or not {rjx^r = 0). In 
the latter case Na{x) is left unmodified (as far as reaction r is concerned), whereas in the 
former we subtract the number of reacting particles {val) and add the number of particles 
produced (far )• 

The probability with which rjx^r is equal to 1 determines the rate at which reaction r 
takes place and should therefore be chosen to reflect the physical situation. In the physical 
case, a reaction rate depends both on the rate constant {K) and on the concentration 
{Co) of reactants; thus for the reaction a + 26 — > c, the standard rate law gives the rate 
of production of c as the product KCaC^- Therefore, the probability of reaction has to 
depend on the rate constant and on the number of reactant particles available at the site 
of interest. Whereas we treat the rate constant as an intrinsic measurable property of the 
reaction, the functional dependence of the microscopic rule on occupation numbers ought 
to yield the concentration dependence of the standard rate law in the continuum limit. 
Thus, the probability that reaction r will take place at x (i.e. for rjx^r = 1) can be written 
as 

p{Vx,r = 1\{M(,{X) : Sf, G S}) = PrFr{{Mp{x)}), (2.6) 

where S is the set of all species. Pr is a real number in the interval [0, 1] and will be 
related to Kr later. is a function of the occupation numbers and by varying it we can 
define various microscopic rules. 

A simple choice for F^. is a function that ensures the presence of sufficient particles for 
the reaction to take place once at the microscopic level. One occurrence of reaction r means 
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that Va} a-particles are subtracted and added at the particular site. According to 

the rule defined with this F^., if there are sufficient reactant particles present, the reaction 
will take place at most once , with probability Py. Formally we define Er as follows: 

Fr{{Np{x) ■.sp(^S})^X[e (Np{x) - (Rule I), (2.7) 

where the 0-function guarantees that sufficient particles are present and the product runs 
over all species. We recall the convention that v^^ = if S/j is not a reactant in reaction 
r: if that is the case, the contribution of species [5 to the product is trivially equal to 
one, since all AZ/j's are non-negative, so that depends effectively only on the occupation 
numbers of the species involved as reactants in reaction r. 



According to the rule defined by Eq. (2.7), once there are sufficient reactant particles 



for a reaction to take place, the latter proceeds with a probability which does not further 
depend on the occupation numbers. Rule I leads to the standard rate law in the continuum 
limit, but only for low particle densities. On the other hand we expect that, the greater the 
number of particles present, the likelier should be the reaction. We introduce therefore 
another rule, in which the reaction probability is weighted by a product involving the 
number of particles of each species. We allow, namely, species (3 to contribute a factor 
Y{m=i{J^I3{x) — m + 1) to Fr. The principal motivation for this choice is that it leads to 

the standard rate law for any particle density. Thus, the new rule amounts to defining 

(») 

Fr{{Mp{x):sp£S})=\{\{{Mp{x)-m + l) (Rule II). (2.8) 

/3 m=l 

Here we allow the product to run over all species, with the convention nm=i(-^(^) — ^ + 
1) = 1, which ensures that species (3 will contribute a factor 1 if it is not a reactant in 
reaction r {v'^^^ = 0). 

Eqs. (13), dU), (U), (Ul) and or (U) define the rule of evolution of the system 
of interacting particles. In the following subsections we are going to use these equations 
in order to derive discrete evolution equations, analogous to the finite difference equations 
(FDE's) approximating the differential reaction-transport equations. For this purpose we 
need to define a real-valued field, the particle density pa{x,t), from the integer-valued 
occupation number Na{x,t). It is clear of course that only Na{x,t) is involved in the 
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simulation, whereas Pa{x, t) can be optionally defined at the level of the output and is not 
fed back into the simulation. 

The occupation numbers obtained in a simulation of the discrete model possess a gran- 
ularity which is not manifest in a typical measurement. Measured quantities are effectively 
averaged over a volume depending on the spatial resolution of the measuring apparatus. 
It is anyway clear that physical attributes of extended systems can only be defined after 
averaging over a certain minimal length scale. The continuum approach appears to de- 
scribe physical quantities over arbitrarily short distances, but this becomes meaningless 
below the above minimal scale. Our approach lies closer to reality in that the transition 
from integers (occupation numbers) to real numbers (particle densities i) is effected via an 
averaging procedure. This can be performed either over a neighbourhood of lattice sites 
or over corresponding sites in an ensemble of states obtained by several independent simu- 
lations of the system. We make a smoothness assumption, according to which the particle 
density does not vary appreciably on the length scale of the neighbourhood used in the 
first of the above averaging procedures. The average densities of two such neighbourhoods 
may of course differ appreciably if their separation is on a larger scale. Under the smooth- 
ness assumption, the results of the above two averaging procedures should agree, within 
fluctuations. Based on this discussion, we define the particle density 

Pa{x,t)=<Na{x,t)>, (2.9) 

where < . . . > denotes an ensemble average. 
2.2 Transport Equation 

Prom a numerical point of view, transport differs from chemical reactions in two impor- 
tant respects: (i) The former involves communication of information between neighbouring 
sites; in the continuum limit this gives rise to derivative terms (advection, diffusion), which 
depend on the precise way the limit is taken, (ii) In our model, the evolution equations 
for pure transport are linear in the density and the respective FDE's are subject to well- 
defined stability criteria; this distinguishes them from the reaction-transport equations, 

■^We reserve the concept of concentrations for the physically measurable quantities, which will be ob- 
tained later by multiplying the densities by a scaling factor. 
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which are usuahy non-hnear and fohow no general stabihty criteria. For these reasons it 
will be worthwhile to concentrate first on pure transport. We begin, therefore, by replac- 
ing C in Eq. by the identity operator. Then the species propagate independently of 
each other and we can drop, for the purposes of the present derivation, the index a. 
Combining Eqs. (2.9) and (2.1), the latter with C = 1, we obtain 



p{x, t + T)=< N{x, t + r) >=< TN{x, t) > . (2.10) 
According to Eq. (|2l^), 

oo 

<TN{x,t)> = Y,(^<^l-_ln><9{N{x-X,t)-n)> + <C(^l><9{N{x,t)-n)> 

n=l 

+ <d+'A,n>< W^ + A,t)-n) >) 
= p< N{x - \,t) > +{1 - p - q) < N{x, t) > +q< N{x + \,t) > 
= pp{x-X,t) + {l-p-q)p{x,t)+qp{x + X,t). (2.11) 



Eq. (2.11) is exact and follows from the fact that the ^'s are statistically independent of 
the A^'s. In deriving ( |2.1l| ) we made use of the expectation values of the components ^i"Jn: 

f < d~n^ > \ /l\ /0\ /0\ / p 



\ < > J 



--<t.,n>=P\ \+{l-p-q)\ 1 \+q\ 



1 — p 



V 



(2.12) 



as well as of the simple result 

oo oo 

^ < 0(iV(x,t) -n) > = <"^e{N{x,t) -n) > 



n=l 



n=l 



< N{x,t) > 



,Vx G £. 



(2.13) 



The last equality follows from the fact that the ^-function contributes 1 to the sum for 
each value of n between 1 and N{x, t) and for higher values. We shall show the same 
result in a different way later (Eq. ( |2.29| )). 

From Eqs. ( |2.10| ) and ( p.ll| ) we deduce the evolution equation 



p{x, t + r) = pp{x - X,t) + {1 - p - q)p{x, t) + qp{x + A, t). 



(2.14) 



Eq. ( |2.14| ) can be readily rearranged as follows 

p(x,t + T) - p{x,t) 
T 

y Pjx + X,t)- p{x - \,t) ^ ^ p{x + A, t) - 2p{x, t) + p{x - A, t) ^2^^^ 



2A 



A2 
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where we have defined 



V^{p-q)- , D^{p + q)^ . (2.16) 



In Eq. (|2l5|) we recognise the forward difference approximation for the first time deriva- 



tive and the central difference approximations for the first and second space derivatives of 
the density: 

p{x,t + t) - p{x,t) _ dp{x,t) 



T 



dt 



+ 0(t) 



p{x + X,t) - p{x - X,t) _ dp{x,t) 2n 

p{x + X,t)-2p{x,t)+p{x-\,t) _ d^p{x,t) , „.,2^ 
- —d^^^^^'- 

Substituting these approximations in Eq. ( 2.15| ) we arrive at the equation 



^ + 0(r) = + D^4^ + O(A^). (2.18) 

ot ox ox^ 

Taking the hmit A ^ 0, r — > and p — (7 ^ in such a way that A^/r and {p — q)X/T 
remain finite, we obtain the transport equation 

dp{x,t) _ y dp{x,t) ^ ^ d'^p{x,t) 
dt dx dx'^ 

On the RHS of Eq. (|1|) there is an advective term, with velocity V, and a diffu- 



sive/dispersive term, with diffusion coefficient D. We note that the continuum limit 
is taken in such a way that V and D, as defined in Eq. ( 2.16| ), remain finite. Eq. (|]l|) is 



the forward-time centred-space finite difference approximation to the transport equation 
( 2.ig| ). We have thus shown that, as far as pure transport is concerned, our model con- 



stitutes a stochastic way of solving the finite difference approximation to the continuum 
PDE. 

A few remarks are in order here: 

(a) The differential equation obtained in the continuum limit depends on the way the 
limit is approached. Let us first try to keep p — q finite. We consider the migration of 
an ensemble of particles concentrated initially at the same lattice site. The density 
of particles evolves according to Eq. ( 2.15| ) , whose RHS contains a term proportional 



to p — q ('advection') and one proportional to p + q ('diffusion'). After a finite time 
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t ^ T, diffusion would result in an average displacement ~ ^/JFIT/t), whereas the 
particles would propagate a distance ~ Xit/r) in the same time due to advection. As 
A ^ and r — > 0, but keeping A/r constant, an infinite number of updates (given by 
t/r) would be needed to make the latter distance finite, but it is after an infinitely 
larger number of steps (given by (t/r)'^) that the former displacement would have 
a chance to become finite. In other words, diffusion would become infinitely slower 
than advection and the ensemble of particles would move without spreading. Thus, 
if we keep p — q and A/r finite, then only the advective term survives as A, r — > 0. If 
on the other hand one tries to keep p—q and A^/r finite, advection becomes infinitely 
fast and this particular limit is of no practical interest. The two processes will be 
comparable in the continuum limit, however, if p — q vanishes like A. p — q is the 
bias in right /left displacement and gives the size of the advective velocity in units of 
lattice spacings per time step. By letting p — q — > we curb the uncontrolled growth 
of the advective displacement, while preserving a finite average diffusive spread. 

(b) From the definitions of V and D and the obvious conditions p — q < p + q < 1, we 
deduce 

< 2D^ < 1 . (2.20) 
A A^ 



In Eq. (|]20D we recognise the stability conditions for the FDE ( |2.15[ ) [p!^ , Eq. (5.1.18)] 



These stability conditions are imposed by hand in the usual FDE approach to en- 
sure that the round-off error introduced by numerical computation does not increase 
exponentially. Since we have shown the equivalence of our method to the FDE, it 
is not surprising that these inequalities hold, but it is an important feature of the 
random walk that they are implemented in an automatic and natural fashion. 



(c) If p + q = 1, the system is subject to the so-called 'checkerboard parity' ||l^. Thus, 
if we colour sites in a checkerboard fashion, two particles occupying at one time 
sites of different colours will always be on differently coloured sites and will never 
meet. The system divides into two subsystems which alternate between sublattices 
of different colour, but remain forever decoupled. We shall show later in this sec- 
tion that, when particles are placed on a lattice according to a uniform random 
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distribution, the occupation numbers obey a Poisson distribution (Eq. ( p. 30 )). The 



existence of two decoupled sublattices does not influence this result, as long as the 
particles are initially distributed evenly between the two sublattices. The opposite 
extreme would be to place all particles initially on sites of one sublattice and then 
perform a large number of diffusion steps. In that case, one of the two sublattices is 
alternately empty and the occupation numbers end up satisfying the Poisson distri- 
bution as if all particles were uniformly distributed only on the occupied sublattice: 
in Eq. (2.30) the density has to be doubled (or, equivalently, evaluated by dividing 



through half the number of lattice sites) and the probability for occupation number 
has to be augmented by a probability of 1/2 that the site belongs to the empty 
sublattice. In intermediate cases, in which particles are placed with a bias favouring 
one sublattice, there will be a corresponding deviation from the Poisson distribution, 
if the latter is calculated on the assumption that particles are distributed over the 
entire lattice. The symmetry just described will obviously not hold if p + q < 1, 
i.e. if particles have a non-zero probability to remain stationary and thus populate 
the same sublattice on successive time steps. Checkerboard parity can also break 
down because of boundary conditions. Thus, the subsystems may mix if we impose 



periodic boundary conditions |12] 



(d) Species with different transport properties (advection velocity or diffusion coefficient) 
can be described on the same lattice by (i) moving particles belonging to different 
species by different numbers of lattice spacings at each update, or (ii) by moving 
species at different multiples of a time step. This way we can simulate, for instance, 
a problem of pure diffusion, where the diffusion coefficients of different species are in 
the ratio of integers. A more natural and physically appealing way to describe species 
with transport coefficients whose ratio may vary continuously is to make p and q 
species dependent. Given Va and Da, we evaluate pa and qa from the equations 

T VaT T VaT 

Pa = + — - , qa = - — - . (2.21) 

If we choose r and A so that (i) pa + qa = ^ for the species with largest Da ^ and (ii) 



'in principle pa + qa < 1 for the species with largest Da will do as well. This choice, however, slows 
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\Va \ X/2Da < 1 for the species with largest \Va \ /D^, then the conditions > 0, 
> and Pa + < 1 will be fulfilled for all species. Condition (i) determines r 
as a function of A and condition (ii) makes sure that A is small enough to make the 
first term on the RHS of the equations ( p. 21 ) larger in absolute value than the sec- 



ond. We note in passing that one might think of treating cases with inhomogeneous 
transport parameters by making the probabilities position dependent. In that case 
the evolution equation ( |2.14| ) becomes 

p{x,t + T) =p{x-\)p{x-X,t) + [l-p{x) - q{x)]p{x,t)+q{x + X)p{x + X,t), (2.22) 

where p and q are functions of position and the species label need not be explicitly 
indicated. From Eq. ( p. 22]) we readily derive the transport equation 

= [Vix)p{x, t)] + ^ [D{x)p{x, t)] , (2.23) 

where V{x) and D{x) are related to p{x) and q{x) according to Eq. ( p. 16 ). Note 
that Eq. ( p.23| ) can be written in the more familiar form Q 

dp{x,ty 



by defining W{x) = V{x) — dD{x)/dx. 



D{x)- 



dx 



(2.24) 



(e) If we repeat the above derivation for a rectangular lattice in d dimensions, assuming 
the same coefficient of diffusion in all directions, we arrive at the resuh D = PX^/lrd, 
where P is the probability that a particle moves to another site during a transport 
step {P = p + q in one dimension). P is a measure of particle mobility and is 
inversely proportional to the retention factor IZ that will be introduced in Section ^, 
in connection with sorbing species. 

2.3 Chemical Kinetics 

We now consider the full problem in which an evolution step is completed by the action 
of the chemical reaction operator C on the result of the transport operation. From Eqs. 



down the simulation and is not ne cess ary, unless the value of t has to be small for reasons such as those 



explained at the end of Subsection 2.3 
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IJ) and (|2J) we obtain: 

Na{x,t + T) = CoTNa{x,t) 



R 

TN^ix, i) + E - ^^}) Vx,r (2.25) 



r=l 



and 



Pa{x,t + T) = <Na{x,t + T)> 

R 

= < TNa{x, t) > + E {^ir - Pr < Fr{{TNp{x, t) : G S}) >, 



(2.26) 



r=l 



where we have directly substituted for < rjx^r > the expression derived from Eq. ( |2.6| ). 

At this point we need a specific ansatz for in order to proceed further. Assuming 
the form given in ( ^.7| ) or (|2.8| ), we are faced with the problem of taking the expectation 
value of a product of random variables. The simplest possibility is that these variables are 
mutually independent. This will be true if the occupation numbers of different species are 
mutually independent; in other words, if the number of particles of every species at a site is 
not correlated with the numbers of particles of the other species at that site. Correlations 
can arise as a consequence of interactions among the particles, i.e. collisions or reactions 
of all kinds. In our model, particles do not explicitly collide, but are scattered by a 
random background (we can think of each random scattering of a particle as simulating 
several collisions of a solute molecule with solvent or other solute molecules). Therefore, 
correlations can only arise as a result of chemical reactions. We shall see in Section |^ that 
correlations do occur in our model. We postulate, nevertheless, molecular chaos, i.e. the 
absence of correlations, for the purposes of our derivation and return to this point in the 
next section. Then the average of the product equals the product of the averages of the 
individual terms and we have to evaluate expressions of the type (^6 {t Na{x,t) — Val^^ 
and < Y\'^Li{TNa{x,t) — m + 1) >, for rules I and II respectively. 

We first derive the evolution equations using reaction rule I (Eq. (|2.7| )). By definition, 
9 (TNa{x,t) —VaTj is equal to 1 when there are not less than Var a-particles at site x 
after the transport operation and otherwise. It follows that the ensemble average of the 
0-function gives the fraction of ensemble members for which there are not less than Var 
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a-particles at x. We define this fraction to be the probabihty that there are sufficient 
a-particles at that position: 

p (TiV,(x, t) > ^ (e (TN^ix, t) - i/W) ) . (2.27) 

{Na{x,t) : Sq, G 5} is as much a set of occupation numbers as {TNa{x,t) : E 5} and 
we can temporarily neglect T: 

p{NUx,t) >n) = {eiN^{x,t)-n)) , (2.28) 

where n is any non-negative integer. It is easy to see that 

oo oo oo 

J2 PiNaix,t) >n) = P(^a(x,t) = m) 

n=l n=l m=n 



oo 



mp{Na{x,t) = m) 



m=l 



<No,{x,t)>. (2.29) 



We note that Eq. (|1|) follows from Eq. ( ^I29|) and the definition (|]2|). 



We now need to relate the probability p{Na{x,t) = n) to the density pa{x,t). We 
are able to fulfill this task under the smoothness assumption, which guarantees that our 
definition of the density as an ensemble average ( Eq. ( p.9D ) is equivalent to the alternative 
definition as a local spatial average. We employ, for the purpose of the present argument, 
the latter definition and recall that, according to the discussion above Eq. (p.S|), local 
averages are evaluated over sections of the lattice which are inhabited by essentially ho- 
mogeneous populations of particles. Following Ref. [0, let such a subsystem contain, at 
time t, N a-particles distributed evenly on M lattice sites. If site x belongs to the par- 
ticular section of the lattice, then pa{x,t) ~ N/M and p{Na{x,t) = n) can only depend 
on n and pa{x,t). We wish to calculate the probability p{Na{x,t) = n) that n particles 

( n\ 

(out of A^) are found at x [M,N ^ n). There are I j ways of selecting n particles 
out of N. Once the n particles are placed at site x, the remaining N — n particles can 
distribute themselves in (M — l)^~" ways on the remaining M — 1 sites. Thus there is a 
total of ^ ^ (M — 1)^~" configurations with n particles at x. The desired probability is 
obtained by dividing the number of these configurations by the number of all possible 
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configurations of N particles on M sites: 

/ N \ {M - 1)^-" _ N\ (M - 1)^-" 



p{Na{x,t) = n) 



~ " , — -e 



{N -ny.nl 



for n<^^ M,N 



(2.30) 



which is just the Poisson distribution. Adopting, for convenience, the compact notation 
Pa{x, t + r/2) =< TNa{x, t) >, we write by analogy 



nl 



(2.31) 



No overdue significance should be attached to the fraction in the time argument t + t/2: 
it merely denotes an intermediate situation just before execution of the chemical reaction 
step. 

Substituting the ansatz (2.7) in Eq. ( p. 26 ), applying the factorisation of the expectation 
value that follows from the molecular chaos hypothesis, employing the definition ( p. 27 ) 



and finally assuming an equality for Eq. (2.31), we find 

R 

Pa{x,t + T) = Pc,ix,t + T/2) + Y,{'^W -'^^a})PrY[p{TNp{x,t)>U^''^ 



r=l 
R 



/3 



p„(x,t + r/2)+^(4{)-4v))p,n E P('^^/5(^'*) 



n 



r=l 
R 



p„(x,t + r/2) + P. n E ^M^l^±^e-M-*+^/^), 



/3 n=..« 

/3r 



n\ 



(2.32) 



where 



Pq(x, t + r/2) = PaPa{x - \,t) + {l-pa- qa)Pa{x, t) + qaPa{x + A, t). (2.33) 



Rearranging terms as in Eq. (2.15) we deduce from Eqs. (2.32) and ( |2.33| ) 

Pa{x,t + t) - Pa{x,t) 
T 

_ Pa{x + X,t) - pa{x - X,t) ^ ^ pa{x + A,t) - 2pa{x,t) + Pa{x - X,t) 



2A 



A2 



R oo 

+ E (-ii^ - -S) n E ^ + r/2)e-''''(-*+-/^), (2.34) 



r=l 



/3 „=vW 

fir 



where 



A 



Da = {Pa + ' 



2r 



(2.35) 
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and the rate constant kr is defined as 



Pr 



(Rule I). 



(2.36) 



kr should not be confused with the physical rate constant K^- which will be defined shortly. 

This is as far as we can go with the discrete model. We now wish to derive the 
continuum limit of the full model. With the continuum approximations (2.17), Eq. ( p. 34 ) 
becomes 



dpa{x,t) 

dt 



+ 0{t) 

R 



dx 



+ 



■r=l 



,(/) 



/5 n=4") 

pr 



(2.37) 



In Subsection 2.2 we let A, r and pa 



0, while keeping A^/r and [pa — qa)^/^ 



finite. In the presence of chemical reactions, we also take Pr — > 0, but keep Pr/T finite for 
all r. It follows immediately that pa — Qa and Pr are 0(A) and 0(r) respectively. Taking 
the continuum limit of Eq. (|2.37| ) as described, we derive the reaction-transport equation 



dpa{x,t) 
dt 



dpa{x,t) d^Pa{x,t) 



dx 



+D, 



dx"^ 



R 

+E 

r=l 



ar ar 



/3 



oo W| 



pr 



p^(x,t)e-'''5(^'*) (Rule I). 



(2.38) 



We note that the forward-time centred-space FDE corresponding to Eq. (p.38|) is similar 



but not identical to the discrete evolution equation (2.34). The difference lies in the time 
at which /o^ is evaluated: for the FDE one uses the value of p/^ at time t, while for 
the evolution equation p/^ is evaluated at the end of the transport step performed at t. 



Moreover, iteration of Eq. (2.34) involves a two-step process where pp{x,t + t/2) is first 
generated through Eq. ( |2.33| ). 

To obtain the standard form of the rate terms we have to consider low particle densities 



{Pa ^1)- If the latter are sufficiently low, then the sums in Eq. ( 2.38| ) are dominated by 



■1— r / 

the lowest-order terms and we obtain products of the type Yip Pp 



pr . 



dpa{x,t) dpa{x,t) d'^pa{x,t) 
— — Vn h L>. 



r=l 



Eq. ( p. 39 ) is typical of the kind of equations numerically solved by conventional ap- 
proaches. It is therefore desirable that the discrete model reduce to a set of such equations 
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in some limit, under well-defined assumptions. If that is the case, it is sensible to test 
simulations of the discrete model by making sure that their predictions converge to the 
solutions of the appropriate differential equations. The chemical reaction rule defined by 
Eqs. ( |2.6[ ) and ( p.7D leads to Eq. ( 2.3S| ) only for low particle densities. As we shall 



demonstrate in Section ^ this restriction to low densities severely limits the efficiency of 
simulations that use Rule I. 

We now employ rule II and show that it leads to the standard reaction-transport 
equations independent of particle density. We substitute Eq. ( p.Sp in Eq. ( p. 26 ) and 



calculate the expectation value of Fr{{TNp{x^t) : sp G 5}). If we assume, as before, 
molecular chaos, the expectation value of the product over (3 factorises. If we further invoke 
the smoothness assumption and the Poisson distribution which this implies ( Eq. ( 2.31| ) ), 



then the individual terms of the product have the form 
< \{{TN^{x,t)-m + l)> = ^ n (n-m + l)p(r7V„(x,t) =n) 

m=l n=Om=l 

= y ff (n - m + 1) ^ + ^/^^ p-P.(^.t+r/2) 



-e 

^^ 1 



(i) In — far 



/2) 



V ^/ H) [n-Var)\ 



' 71=0 



= + (2.40) 

On the RHS of the first equality above, the summation over n need only be carried out 
from n = z/ar , because the product over m always vanishes for < n < Var — 1. Repeating 
the steps of our earlier derivation of the continuum limit, we arrive at ( cf. Eq. ( |2.38| ) ) 

^ = -^"^ + + E - S) KRpfi-.') (R"le II), 

r=l f3 

(2.41) 

where 

kr = Pr/r (Rule II). (2.42) 
Eq. ( ^ ) holds as a strict equality without the additional assumption of low density. 
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If Pa were solute concentrations, this would be the form commonly used in modelling 
systems with transport and chemical reactions. We have already remarked that particle 
densities are not the same as measurable concentrations. Clearly, we do not expect in 
the foreseeable future to be able to treat particle numbers comparable to Avogadro's 
number. We therefore assume that it is legitimate to work with numbers of particles that 
are by orders of magnitude smaller than Avogadro's number and obtain the respective 
concentrations Ca{x^t) by rescaling the densities pa{x,t): 

Ca{x,t) =-fPa{x,t). (2.43) 

We further assume that the scaling factor 7 is universal (i.e. independent of species, space, 
time and the value of the density itself). Thus, we typically fix the value of 7 from 
the initial conditions of the simulation (total particle density) and the real problem it 
purports to model (total concentration). We then perform the simulation and recover the 
concentrations at the desired times and locations by multiplying the densities at those 
times and locations by 7. If we substitute pa = Ca/j into Eq. ( p. 41 ), we can absorb 7 in 
the rate constants by defining the new rate constants Kr- 

Kr = kr-f'^c,''^^r+\ (2.44) 

Then, the concentrations satisfy equations of the standard form: 

r=l f} 

Here the parameters Va,Da and model directly properties of the physical system. In 
particular, is identified with the physical rate constants in Eq. ( ^.4] ). We note that, if 
we substitute Pa = Call into Eq. (|2.38 ), then 7 cancels out if only linear terms (transport 



terms and, possibly, linear reactions) are present, but in general non-linear terms result in 
a non-trivial 7-dependence. For low densities, rule I yields Eq. ( ^.45 ) in an approximate 
form, following from Eq. (1^. 



It should be mentioned that even rule II is not completely free of limitations on the 
density. The reason is quite different this time and becomes obvious if one substitutes 
(|2.8| ) in Eq. ( |2.6D . Since the probability for rix^r = 1 rnust be < 1 (we refer to this 
condition as probability conservation), the product of Pr and the appropriate combination 
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of occupation numbers at any site should not exceed one. For example, if r is the reaction 
a + 25 — > c, we should have Pj-NaNfj^Nfy — 1) < 1. One has to make sure, therefore, that 
particle density is low enough to guarantee that the probability for the above product to 
exceed one is negligibly small. Alternatively, since it is the ratio Pt/t that is related to 
physical parameters, we can handle arbitrarily high particle densities by making Pr and, 
consequently, r appropriately small. The limit in this case is set by the resulting increase 
in the number of iterations and the accordingly greater computation time. 

In Section |^ we shall study in detail a reaction-diffusion system in which particles of 
the species a, b and c interact via the reversible reaction 

a + b ^ c. (2.46) 

K2 

The reaction-transport equation for the concentration of species a is obtained as a special 
case of Eq. ( p5|) : 



dCa{x,t) dCa{x,t) d'^Ca{x,t) +\ , n f A (O A7\ 

di ^ ~ dx dx^ KiCa{x, t)Cb[x, t) + K2Cc[x,t). (2.47) 



According to Eqs. ( 2.36| ) and ( 2.44 ), if there is at least one a- and one 6-particle at a site, 
the reaction a-|-6 — > c will occur with probability Pi = t'jKi. Similarly, if there is at least 
one c-particle, it will disintegrate into an a and a b with probability P2 = TK2. If we are in 
a diffusion-limited regime, where chemical equilibrium is attained on a much shorter time 
scale than that of diffusion, then a dynamic equilibrium is established locally between the 
reactions a -|- 6 — > c and c — > a + 6, the two rates cancelling each other in Eq. ( p.47| ): 

KiCaCb = K2C,^^ = ^^^ = ^ , (2.48) 

CaCb K2 PaPb k2 



where we have omitted for simplicity the space and time arguments. Eq. ( |2.48| ) is a special 
case of the law of mass action for ideal solutes, with equilibrium constant fC = Ki/K2 - i 
If we use rule I, the exact reaction-transport equations are obtained as a special case 
of ( p.38|) , for example, 

9_PJ^ = (1 - e--(-*)) (1 - e--(-^))+^, (1 - e--(-*)) 



(2.49) 



^The law of mass action in this simple form holds only for infinite dilution. For higher concentrations, 
interactions among the species (e.g. of electrostatic nature) complicate the situation. Such effects may in 
principle be incorporated in our model, but are not considered in its present form. 
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where ki = 'jKi and k2 = i^2- We note that no low density assumption has been made 



here. According to ( 2.49 ), the law of mass action (Eq. ( 2.48 )) is replaced by 

1 - e-P'^ _ki 
(1 - e-P-) (1 - e-Pb) ~ Y2 ■ 

For Pb, Pc ^ ^ recover the familiar law of mass action and, if we write Ca{x,t) = 
^Pa{x, t) etc., Eq. ( 2.4g| ) reduces to ( 2.47| ) as an approximate equation; if we are sufficiently 



close to chemical equilibrium, KiCaCb = 'ykipaPb is of the same order of magnitude as 
K2CC = ^k2Pc and we can say that Eq. (\2A7\) holds to 0(0^/7^) = 0{pl). 



2.4 Homogeneous System 

In a homogeneous system, the particle density is independent of the spatial variable 
and we can write it as Pait)- We define Pait) as the average number of a-particles per 
site: 

Pa[t) = , (2.51) 

where ^aix, t) is the total number of particles of species a and V is the total number 
of sites. In two dimensions, with sites in the x-direction and Ny sites in the y-direction, 
V = Nx X Ny. For a homogeneous system, this definition is equivalent to the ensemble 
average used in Section 

We repeat the steps that led from Eq. ( |2.25D to Eqs. ( p. 381 ) and ( |2.41| ), but with < . . . > 



understood this time as the average over all sites. The first term on the RHS of Eq. ( p. 25 ) 



becomes < TNa{x,y,t) >= Pa{t) =< Na{x,t) > upon averaging, as if the transport 
operator T had not been applied at all to the occupation number. This is natural, as 
transport does not change the average properties of a homogeneous system. We factorise 



the expectation value of Fj. in Eq. ( |2.26 ) under the assumption of molecular chaos (see 
the discussion following Eq. ( |2.26| )) and evaluate the individual terms using the Poisson 
distribution (derived on the basis of the arguments that led to Eq. ( |2.3C| ) , but applied this 
time to the system as a whole). We thus arrive at an equation for pait). 

^ = E {-ii^ - n E -^Ppit)--'"^'^ (R^le I), (2.52a) 

^ = E (-ii^ - -S) kr n pf (0 (Rule II). (2.52b) 



r=l 
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2.5 Boundary Conditions 

We shall discuss here the transport behaviour of particles when they reach the lattice 
boundary. Since chemical reactions will not enter in the present discussion, we shall make 



use only of Eq. (2.10). Each boundary condition will be determined once we formulate 
the microscopic transport rule which replaces Eq. ( |2.2| ) for sites at the lattice boundary. 

We first consider an impermeable boundary, with particles bouncing back when they 
reach it. In one dimension we can set the boundary at the first site on the left (say, 
X = 0) and the last on the right (x = Xmax)- We consider the left end of the lattice. The 
first two sites are x = and x = X respectively. We define the transport operation as 
before by the rule: $xn = (1)0) 0), (0, 1, 0) or (0, 0, 1) with probability p,l — p — q and 
q respectively, unless x = 0, in which case ^xn = (1)0)0) or (0,1,0) with probability p 
and 1 — p respectively {p and q are the same as in the interior of the lattice). In other 
words, once a particle reaches the left boundary, it moves to the right with the same 
probability p as in the interior of the lattice and remains at the boundary with probability 
1 — p. Equivalently, one may think of the boundary as lying at —A/2 with particles at 
X = moving according to the same rule as in the interior, but bouncing back to x = 



within the same transport step if they hit the boundary. The equivalent of Eq. (2.11) for 
X = is obtained from the above probabilities and the fact that there are no particles 
coming from the left: 

p(0, t + T)=< TN{0, t) >= (1 - p)p{0, t) + qp{X, t), (2.53) 



Eq. ( |2.53| ) can be rearranged as follows 

p(0,i + r)-p(0,t) _ qp{X,t)-pp{0,t) 



T 



A/. , . p{X,t) - p{0,t) p-g . . . fn 
-iip + Q) ^ ^ [p(A)t) + p(0,t)]| 



(2.54) 



In the limit A, r ^ 0, with p - q 0(A) and A^/r ~ 0(1), the LHS of Eq. ( |2.54| ) remains 



finite while on the RHS A/r — > oo. This implies that the expression in the curly brackets 
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on the RHS vanishes so that 

= [P(A, t) + p{0, t)] = ^ [p{X, t) + p{0, t)] , (2.55) 



where we have used the definitions (2.16). In the continuum hmit we obtain 

dp{x,t)' 



Vp{x,t) - D- 



dx J ^=0 



. (2.56) 



Eq. ( 2.56| ) is a statement of the condition that the total flux of solute (i.e. the sum of the 



advective flux Vp{x, t) and the diffusive/dispersive flux —Ddp{x, t) jdx ) vanishes at x = 0. 
This is intuitively clear from the definition of the impermeable boundary: any particles 
that reach it bounce off so that there is no flux across the boundary. 

It is common in solute transport problems to specify either the concentration or its 
gradient at the boundary. These or mixed boundary conditions ( i.e. relating the concen- 



tration with its gradient, such as Eq. ( 2.56| ) ) can be easily implemented in our simulations 



once their physical background is clear. Thus, in the case when the concentration at the 
X = boundary is fixed, it may be assumed that there is a large homogeneous reservoir, 
extending beyond the system of interest and having the given concentration. Introducing 
explicitly a reservoir beyond the boundary at a; = would be correct but impractical. 
This situation can be more efficiently simulated by assigning to x = at each transport 
step an occupation number from a set of random numbers obeying the appropriate Poisson 
distribution (instead of applying the simulation rule at x = 0). A special case is that of 
a sink, i.e. vanishing concentration, at x = 0. The occupation number is then set at all 
times equal to zero at x = and the evolution rule is applied normally to all sites of the 
lattice interior. 

Alternatively, one may specify 5/j(x, t)/5x|^^Q = as the boundary condition. To 
simulate this, we add formally site x = — A to the lattice. Before each transport step we 
set the occupation number at x = —A by N{—X,t) = N{X,t). The evolution rule is then 
applied to all other lattice sites, including x = 0: 

p{0,t + T) = pp{-X,t) + {l-p-q)p{0,t) + qp{X,t) 

= pp{X,t) + {l-p-q)piO,t)+qpiX,t) (2.57) 
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which can be rearranged as 

p{^:t + T)-p{Q,t) ^^^^ ^^ p{X,t) - p{0,t) ^ 

T r 

Arguing as before, we must have p{X, t) — p{0, t) ~ O(A^), if the RHS of ( p. 581 ) is to remain 
finite in the continuum hmit. It follows that dp{x,t)/dx\^^Q = [p{X,t) — p{0,t)] / X + 
0(A) ^ 0, as A ^ 0. In physical terms, we make the diffusive flux vanish at the boundary 
by superposing equal and opposite amounts of outgoing and incoming diffusive flux; this 
leaves only advection to take care of net solute transport across the boundary. 

A periodic boundary condition is often computationally convenient and is used when 
the precise behaviour of the boundary layer is not important, e.g. with translationally 
invariant homogeneous systems or when the boundary is too far away to influence the 
region of interest. We impose a periodic boundary condition by connecting the two ends 
of the lattice, so that particles crossing the right boundary appear automatically on the 
left boundary and particles crossing the left boundary appear on the right boundary. 
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3 Simulation of reaction-transport processes 

The model formulated in Section |^ is very general and can describe a wide variety 
of coupled transport-chemical reaction processes. In this section our primary aim is to 
demonstrate the versatility of the model by using it to simulate various systems. After 
demonstrating a simple system undergoing both diffusion and advection, we concentrate 
on systems of diffusing particles subject to various chemical reaction schemes. A thorough 
discussion is given of the reactions a + b ^ c and a + b ^ c; these reactions are particularly 
suitable for displaying the essential microscopic aspects of the simulation. In each case, the 
results of the simulations are compared with the ones obtained from the corresponding 
PDE's i. The purpose of the comparison is to test the validity of the smoothness and 
molecular chaos assumptions used to derive the PDE's in Section |2|, as well as to study 
how the simulation converges to the continuum result. Lastly, we show that our model is 
also capable of simulating complex autocatalytic reaction-diffusion systems which, under 
non-equilibrium conditions, display remarkable spatial and temporal structures. 

3.1 Diffusion and Advection 

We first discuss solute transport without chemical reactions. Advection and diffusion 
arise as the macroscopic result of a random walk. In fact, from a microscopic point 
of view, there is no fundamental difference between the two processes. If the random 
walk is unbiased (equal probability of motion in all directions) there is only diffusion and 
no advection, whereas a bias in favour of a certain direction produces diffusion coupled 
with advection in the chosen direction. Fig. |l| shows the result of a simulation of solute 
transport on a one-dimensional lattice. We simulate simultaneously two solutes with 
advection velocities Vi = 0.5 m/y and V2 = Vi/W and diffusion coefficients Di = 2bm? /y 
and D2 = -Di/10. Taking r = 10~^ y and pi + = 1, we use the second of Eqs. ( p.l6[) 

to determine A. The probabilities pi, qi, as well as ^2, 0.2, are calculated from Eq. ( p. 21 ). 

^The differential equations are solved using standard finite difference methods. We do not intend to 
discuss here the convergence of the standard methods. All results obtained with them have been checked for 
convergence and are indistinguishable for our purposes from the true solutions of the continuum equations. 

^The second solute can be thought of as being retarded by a factor TZ — IQ due to sorption on the 
surface of the solid matrix. If one assumes instantaneous sorption equilibrium and a linear relationship 
between liquid and solid phase concentrations, the effect of sorption can indeed be reduced to a retention 
factor TZ that divides the transport coefficients Va and Da- 
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Figure 1: Result of one-dimensional simulation (species 1: solid circles, species 2: open 
circles) compared with the solution of the transport equation (species 1: solid curve, species 
2: dashed curve) at t = 240y. Parameters and initial conditions are given in the text. 

Initially both concentrations are equal to 1 (arbitrary units) for x < 0, and for x > 0. 
The simulation begins with 200 particles of each species per site to the left of x = and the 
concentration is calculated by averaging the occupation number over cells of 100 lattice 
sites and normalising by a factor 7 = 1/200. The boundary condition is uninteresting 
here, because the boundary is chosen far enough, so that its influence does not reach 
the displayed region by the time considered. In the figure, the solution of the transport 
equation ( 2.19] ) is shown for comparison. Error bars of length equal to one standard 



deviation were estimated and were always found to be smaller than the plotting symbol. 
The small fluctuations around the solid curve are of statistical nature and diminish if we 
increase the number of particles used in the simulation. 

3.2 a-|- b # c , homogeneous case 

We introduce chemical reactions by looking first at a system of particles which is ini- 
tially homogeneous. We begin the simulation by placing particles on a two-dimensional 
lattice according to a uniform random distribution. As we saw in Section ^, for a ho- 
mogeneous system the spatial derivatives of the density vanish and the reaction-transport 
equations ( 2.38 ) and ( 2.41]) , corresponding to reaction rules I and II respectively, reduce to 
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the chemical rate equations ( |2.52a| ) and (2.52b), which describe the evohition of uniform 



particle densities. It is important to note that these equations have been derived under 
the assumption of molecular chaos. In reality, correlations between particles and density 
fluctuations do occur, and on occasion, notably in some irreversible reactions as well as 
in some autocatalytic systems maintained far from equilibrium, lead to inhomogeneities, 
even for systems which initially are homogeneous. Later we shall discuss such reaction 
schemes; then of course, the average particle density need not follow Eqs. ( |2.52a| ) and 
(2.52b). For the moment, however, we consider the case where the system remains homo- 



geneous as a function of time, and we address the question of how microscopic dynamics 
drives a discrete system towards equilibrium. 

(a) Reaction rule I 

We consider a homogeneous system of a-, b- and c-particles reacting via the reversible 
reaction a + b ^ c. At each time step every particle moves randomly to one of the four 
nearest sites. Particles react according to rule I: if there are at least one a- and one b- 
particle at a site, then, with probability Pi, one a- and one 6-particle are removed and 
a c-particle is added; whereas, if there are one or more c-particles at a site, then one of 
them is replaced by an a- and a 6-particle with probability P2. As long as we discuss 
homogeneous systems, we shall average densities over the whole lattice and express them 



only as functions of time, Pa{t) (cf. Subsection 2A). We follow the approach to equilibrium 
by looking at the time evolution of the reaction quotient pc{t) / pa(t) pb{t) . 
Let 

Pax {Qax) be the probability that an a-particle moves by one lattice spacing A to 
the right (left) along the x-direction and p^y (qay) the probability that it moves up (down) 
by the same distance along the y-direction in one transport step. Assuming all species to 
have the same diffusion coefhcients, we take pax + Qax + Pay + Qay = Va, i.e. particles 
always move to a neighbouring site. 

Here, as in all cases below, we consider the case of no advection for convenience; when 



necessary, advection can be easily included in our model as demonstrated in Subsection 3A 
above. Thus we set Vax = {Pax—Qax)^/T = and Vay = {Pay—Qay)^/'T' = 0, where Vax and 
Vay are the advection velocities in the x- and y-directions respectively (thus pax = Qax and 
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Pay = Qay, Vo). We further assume equal diffusion coefficients in the x- and y-directions 
{Dax = {Pax + 9a2:)A^/2r = Day = {pay + Qay)^^ /'^t) ■ Putting together the constraints 
on the various displacement probabilities, we deduce pax = Qax = Pay = lay = 1/4 
and Dax = Day = A^/4r = D. For convenience, these conditions will hold for all two- 
dimensional systems in this paper. We begin the simulation with 50 000 particles of each 
species on a two-dimensional lattice of 500 x 500 sites (i.e. a density of 0.2 particles of 
each species per site) with periodic boundary conditions. The reaction probabilities are 
determined through the relations kr = Pr/T, where the rate constants are taken to be 
ki = 0.8 and k2 = 0.2. 

For our initial simulation we take r = 1 (arbitrary units) and examine the time develop- 



ment of the reaction-diffusion system up to time t = 20. In Fig. ^ the solid circles denote 
the mean value of the reaction quotient for an ensemble of 21 systems. The estimated 
error bars are smaller than the plotting symbol. The solid curve is obtained by solving the 
rate equations ( 2.52a| ). According to the latter, the reaction quotient reaches the equilib- 
rium value 3.50 at roughly t = 10 and remains constant thereafter. This value can also 
be derived by solving Eq. ( 2.5C| ) subject to the constraints Pa{t) — Pb{t) = constant and 
Pa{t) + Pc{t) = constant, which obviously hold for a + b ^ c. The result of the simulations 
evolves instead towards a slightly higher equilibrium value and then runs parallel to the 
solid curve. The simulations lead to an equilibrium state with relatively more c-particles 
than predicted by the rate equations. 

That the discrete result does not agree with the exact solution of the rate equations is 
due to two reasons. As far as the rate equation is concerned, the problem is completely 
specified by the initial conditions, the rate constants kr, and the final time of interest, 
t = 20. Apparently, the time discretisation used in the simulation is too coarse to pro- 
vide agreement with these rate equations. To illustrate this, we solve the coupled FDE's 
obtained from Eqs. ( p.52a| ) with the same time step, r = 1, as used in the simulations & 



The result of the FDE's (dashed line in Fig. 2a) lies far from the exact result it is approx- 
imating (showing that the time discretisation is not sufficient); but it also lies far from 
the result of the simulations (although the latter use the same time discretisation as the 



The FDE's are obtained by substituting [pa(i + r) — Pa{t)] /t for dpa(t)/dt in (2.52a). 
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Figure 2a: Time dependence of the reaction quotient for a two-dimensional, homogeneous 
system of a-, b- and c-particles, reacting via a + b ^ c according to rule I: mean value 
over 21 simulations with r = 1 (solid circles), exact result of differential rate equations 
(solid curve) and approximate result of FDE's (dashed curve). Also shown is the result 
of 21 simulations with r = 1/5 (open circles). The values of the transport and reaction 
parameters, as well as the initial conditions, are given in the text. 



Figure 2b: The same as Fig. p^, but with 10 d.p.r. 
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FDE's). 

The discrepancy between the simulations and the solution of the FDE's is of a different 
origin and can be understood as follows: The rate equations and the approximating FDE's 
rely on the assumption that a- and 6-particles are uncorrelated. In the simulations, this 
assumption does not hold. At the rate of one diffusion step per reaction step (d.p.r.), 
the a- and 6-particles originating from the disintegration of a c-particle remain sufficiently 
close to each other for the probability of their meeting again (with a subsequent chance 
of reaction) to be significantly higher than in a uniform distribution of particles, i.e. if the 
assumption of molecular chaos were valid. The reverse reaction, c ^ a + b, is not affected 
by correlations. As a result, a + 6 — > c is enhanced with respect to c — > a + 6 and relatively 
more c-particles are produced, resulting in a higher reaction quotient. The conditions of 
molecular chaos can be systematically approached by performing an increasing number, 
n£), of d.p.r. In each evolution step of our simulation, we repeat the transport operation 
riD times before we perform the reaction operation once. This multiplies the diffusion 
coefficient by a factor n^) in the continuum limit, but, for a homogeneous system, the rate 
equations remain the same. The reaction quotient in the asymptotic state is shown in 
Fig. ^ for riD = 1, 2, 3, 5 and 10. To obtain the value for ud = oo, we replace the transport 
operation by a random redistribution of particles §. We see that, as increased diffusion 
washes out correlations between the particles, the equilibrium reaction quotient approaches 
the value obtained from the rate equations. Similar conclusions concerning correlation 
effects have been drawn in Ref. [|l^]. 

The time evolution of the reaction quotient for ud = 10 is shown in Fig. ^ (solid 
circles). As expected from the preceding paragraph, correlations play a much lesser role 
and the system reaches a steady state only slightly higher than that predicted by the rate 
equations. However, for times t < 5, there is a serious discrepancy with the exact solution 
of the rate equations. Comparing again with the solution of the FDE's (dashed curve), we 
find good agreement, thus confirming the smallness of residual correlations. In fact, the 

agreement becomes perfect if we redistribute particles randomly between reaction steps 

^We have tested the validity of this identification by letting a system of particles diffuse from different 
initial distributions. After sufficient iterations of our diffusion algorithm, the occupation number obeys 
essentially the Poisson distribution over the accessible lattice sites, as expected from a uniform random 
distribution of particles. 
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Figure 3: Reaction quotient for the system of Fig. averaged over 9000 time steps in 
the steady state of a simulation, for different numbers of diffusion steps per reaction step. 



{UD = oo). 

We now look at the way our discrete model approaches the continuum limit. As we 
refine the time step, we expect to obtain the exact prediction of the rate equations, in 
very much the way the solution of an FDE converges to the exact solution. We perform 
simulations with a shorter time step r = 1/5. Then we also have to scale down the reaction 
probabilities Pi and P2 by a factor 5, in order to obtain the same rate constants ki and k2- 
The results of the simulations executed with these parameters are shown as open circles 
in Fig. = 1) and in Fig. ^ (n/j = 10). The reaction quotient is here averaged over 

21 independent simulations and statistical errors are smaller than the plotting symbol. 

The finer discretisation obviously leads to better agreement with the exact result in 



both Figs. 2a and 2t. In the case of Fig. 2b, where correlations are almost absent, the 
improvement over the previous set of simulations (solid circles) is attributable to the use 
of a smaller time step, in the sense of convergent FDE schemes. In the case of Fig. pa| , 
however, the finer discretisation improves the agreement also partly by reducing indirectly 
the effect of correlations. This can be seen as follows: As the reaction probabilities Pr 
decrease (by the same factor, say m, as the time step), particles react more weakly in a 
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single update. This is compensated by exactly that many more updates within a given time 
span, so that the overall reaction strength (rate constant) remains intact. The particles 
execute, however, m times more movements to neighbouring sites in a given time span. 
As a result of this increased agility, particles originating at the same site are more inclined 
to lose memory of their common origin and correlations are subsequently weakened. In 
fact, as the continuum limit is approached, any correlations between the a- and b-particles 
will disappear. 

Using reaction rule I (Eq. ( p.52a| )), we obtain a reaction quotient of 3.50 at chemical 



equilibrium, for the particular parameters chosen. The standard rate law, Eq. ( 2.52b ), 
leads to a significantly different value of ki/k2 = 4.00. This difference is obtained for an 
initial particle density of 0.2. It is clear that the particle density has to be significantly 
lower than 0.2, if we intend to reproduce the standard law by using Rule I. In fact a 
density of 0.005 would be necessary to reduce the discrepancy at equilibrium below 1%. 
For low particle densities the statistics is poor for each individual simulation and one 
has to perform many of them or use a bigger lattice. On the other hand, one can also 
obtain better statistics if higher particle densities can be used. We therefore proceed to 
investigate reaction rule II, which allows us to use higher densities. 

(b) Reaction rule II 

We first apply rule II to a system of particles with the same initial density, 0.2, as with 
rule I. Fig. ^ shows the result of 21 simulations, on a 500 x 500 lattice, of a homogeneous 
system of a-, b- and c-particles diffusing with = 1 and reacting via a + b ^ c, according 
to rule II. Transport and kinetic parameters, as well as initial and boundary conditions, 
are identical with those of the previous simulations on such a lattice. The time step is 
taken to be r = 1/10 and the reaction parameters Pi = 0.08 and P2 = 0.02. According to 
rule II, if there are Ua a-particles and Uf, fe-particles at a site, then one a- and one 6-particle 
are removed and a c-particle is added with probability PiUaUf,, whereas, if there are Uc 
c-particles, one of them is replaced by an a- and a 6-particle with probability P2nc- A 
smaller time step is dictated by the need to make Pi and P2 sufficiently small so that the 
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Figure 4a: Time dependence of the reaction quotient for a two-dimensional, homogeneous 
system of a-, b- and c-particles, reacting via a + b ^ c according to rule II: mean value over 
21 simulations (open circles) and exact result of differential rate equations (solid curve). 



Figure 4b: The same as Fig. 4a, but with an initial density of 1 and appropriate rescaling 



of parameters and densities as explained in the text. 
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probability conservation conditions ( cf. Eqs. ( |2.6[ ) and (|2.8| ) ) 

PiVx,! = l\{Na{x,t) = Ua}) = PlUaTl}, < 1 (3.1a) 

and p{'nx,2 = M{Na{x,t) = na}) = P2nc<l- (3.1b) 

will be respected in all but an insignificant number of cases. For Pi = 0.08, condition 
( 3.1a| ) is violated if Uarih > 12. Assuming a Poisson distribution for the occupation num- 
bers, we readily estimate that naUh > 12 occurs at an arbitrary site with probability 
0.8 X 10"^ II . From this we also estimate that there is a 0.2% likelihood that probability 
conservation will be violated at all (i.e. anywhere on the lattice) in one iteration step. This 
likelihood is negligible for practical simulations and, in the rare cases when the probability 
of reaction does exceed one, we can safely set this probability equal to one. Due to the 
substantially smaller time spacing (r = 1/10), correlation effects are insignificant even for 
1 d.p.r. in this case. 

Moving to a higher particle density, we simulate, on a 200 x 200 lattice, a homogeneous 
system of a-, b- and c-particles, reacting via a + b ^ c, with an initial density of 1 particle 
of each species per site and periodic boundary conditions. The time step is r = 1/10, as in 
the last set of simulations, and we perform 1 d.p.r. To obtain the same rate equations and 
initial conditions as before (and hence the same continuum result to compare with), we 
have to first rescale the previous value of Pi by a factor 7 = 0.2 (leaving P2 intact), then 
perform the simulation and finally multiply the resulting densities by the same factor. The 
scaling factor 7 is a dimensionless version of the factor used to relate particle density to 
concentration in Section ^. The result of 21 simulations is shown in Fig. |4^. In this case, 
probability conservation is violated if UaUb > 62, which occurs with probability 0.3 x 10^^ 
at a given site and with probability 0.1 x 10~^ anywhere on the lattice during an iteration 
step. 

(c) Optimisation of Rule II 

We have seen that rule I approximates the standard rate law only for sufficiently low 
particle densities. Under this constraint, the quality of statistics can be improved either 
by performing more simulations, or by increasing the number of lattice sites. Rule II, 
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however, leads to the standard rate law without the restriction of low particle densities 
and, in that case, better statistics can also be achieved by utilising the higher densities. 
On the other hand, rule I is simpler and its simulation computationally more efficient 
for a given density. It is very important to assess whether the additional freedom (wider 
density range) afforded by rule II can be exploited in order to obtain higher computational 
efficiency (for the same statistics). More precisely, we address the questions: (i) whether 
there is an optimal density which minimises the computation time needed to perform a 
simulation using rule II, (ii) how this computation time compares with the time required 
by rule I. We shall not attempt to provide a universal answer to these questions; instead, 
we concentrate on the example of a + 5 ^ c we have been considering so far, emphasising 
the conclusions that have wider validity. 

We perform a series of simulations of the above system of a-, b- and c-particles on a 
200 X 200 lattice, using rule II and different values of the initial particle density, which we 
take to be the same for all species: Pa{0) = Pb{0) = Pc{0) = P- The reaction parameters 
Pr are rescaled each time so that the continuum rate equations and initial conditions are 
the same as those in the last set of simulations above. We finally rescale the densities 
Pa{t) obtained from the simulation by the ratio of the reference density 0.2 to the initial 
density of particles on the lattice, p, thus Pait) — > Pa{t) = {0.2/ p)pa{t). In each case 
we calculate the reaction quotient pc{t) / pa{t)ph{t) and compare its asymptotic value with 
the result obtained by solving Eq. (|2.52bD with ki = 0.8, ^2 = 0.2, p = 0.2 and the 
stoichiometric coefficients appropriate for a + h ^ c. 

We have already seen that the result of the simulations lies arbitrarily close to that 
of the rate equations for the right choice of r. Here we require that the result of the 
simulations does not deviate by more than ~ 10 — 15% from the exact result. This 
is roughly the discrepancy we obtain if we solve Eq. ( |2.52a| ) with the same parameters 
(resulting in a reaction quotient of 3.50 instead of 4.00, as given by the standard rate law). 
This discrepancy reflects the systematic error expected from rule I, since the densities 
calculated with that rule converge to values that satisfy Eq. ( 2.50| ) instead of the standard 
Eq. ( p^ ). 

Given the above maximum acceptable error, we seek to minimise the computation 
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Figure 5: (a) Normalised minimum computation time for one simulation, as a function of 
initial particle density: Rule II (crosses, curve to guide the eye) and rule I (solid circle), (b) 
Estimated minimum computation time (Tc = c/pPi), corresponding to a 15% likelihood of 
probability-conservation violation at a site, as a function of particle density, with c fitted 
to actual Tc at p = 10. 

time by reducing the number of iteration steps. This is equivalent to increasing r, which 
in turn implies a larger Pi (of course, also a larger P2, but this has no influence on 
most of our considerations). Increasing r can have a number of possible effects on the 
results of the simulation. Firstly, we expect deviation from the rate equation in the sense 
that the corresponding FDE deviates from the rate equation. However, it can be easily 
checked that the FDE gives the same reaction quotient at equilibrium as the rate equation, 
independently of the value of r. Since we have chosen to use only this asymptotic ratio for 
comparison of simulations, increasing r will not give rise to such numerical convergence 
problems in our case. On the other hand, an increase in the probability of reactions leads 
to stronger correlations between the particles, and this does affect the reaction quotient at 
equilibrium. By varying the number of diffusion steps per reaction step, we have actually 
found that the effect of such correlations is relatively small. We find instead that the 
main effect of increasing r, and hence Pi, is that probability conservation is violated more 
frequently. Since an unphysical probability (> 1) is treated in our simulation as if it were 
1, the reaction a + b ^ c is partly suppressed each time Pi exceeds 1. This amounts to 
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a suppression of the production of c-particles and appears macroscopically as a reduction 
of the reaction quotient. 

For a given initial density p, we obtain our minimum computation time, tc{p), when 
the reaction quotient is reduced by the maximum amount allowed, namely by the factor 
10 — 15% discussed above. For all densities we normalise this minimum computation 
time by demanding the same statistics as for p = 0.2; thus we introduce the normalised 
minimum computation time Tc{p) by the equation Tc{p) = jtcip), where 7 = 0.2/p. In 
Fig. |5|a we show Tc as a function of the initial density p: there is a minimum for p between 3 
and 5. Thus, there is a range of densities for which rule II can be simulated with maximum 
efficiency. To compare this with the efficiency of rule I, we also show the computation 
time for a simulation of that rule, performed with initial density 0.2 and the same number 
of iteration steps as the optimal simulation of rule II with the same density. We see that, 
although rule I may be significantly more efficient than rule II for the same density, it can 
be less efficient than rule II when the latter is used with a higher density. In other words, 
using rule II with higher densities can offer an advantage in computational efficiency (in 
the particular case we are considering, by a factor of 4-5). 

In Fig. ^ we do not continue to higher densities because probability conservation 
begins to be violated due to the size of P2 and the picture becomes more complex beyond 
p ~ 10. We expect however that the upward trend of the curve in Fig. ^a continues 
to higher densities. To show that higher densities result in ever increasing computation 
times we argue as follows: It is plausible to assume that, for high particle densities p, the 
computation time Tc required for given statistics is roughly proportional to Ng x Nt x Np, 
where Np is the total number of particles, Nt the number of iteration steps in a single 
simulation and Ng the number of independent simulations required to achieve the desired 
statistics. Since Ns oc 1/p and Np ~ pN^ (assuming a rectangular Nx x Ny lattice, with 
Nx of the same order of magnitude as Ny), we have Tc oc N^ x Nt. But Nx = L/X and 
Nt = T/t, where L is the linear size of the system and T the time interval during which 



the system evolves. For a reaction of the type a + 6 — > c, we saw below Eq. (2.47) that 



r = Pi/^Ki, where 7 relates p to the physical concentration C: C = ^p (cf. Eq. ( |2.43| )); 
from Eq. ( |2.16|) we also find A oc ^Dt. Putting things together, we deduce Tc oc l/(pPi)^. 
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Since our aim is to minimise computation time, Pi has to be as high as possible, with 
the upper bound set by the maximum acceptable likelihood of probability-conservation 
violation. Given the density, we can calculate the violation probability for any given value 
of Pi assuming a Poisson distribution of the occupation numbers. Conversely, if we fix the 
violation probability, we can determine the corresponding value of Pi. To be definite, let 
us set the violation probability at a site to be 15%. It is then found that Pi decreases, as a 
function of p, faster than 1/p, which implies that and hence T^, is an increasing 

function of p. The decrease of Pi, and hence of r and A, with p results in an increase 
of lattice size and iteration number which in turn implies an increase in computation 
time. There is a subtle point in the above argument which deserves mentioning. For a 
homogeneous system we need not require a fixed diffusion coefficient D and hence we can 
keep A and fixed while we vary Pi and r. This is in fact what we did in the the 
simulations described above. Then Tc oc Nt oc 1/pPi, which is shown in Fig. 5b, is also an 
increasing function of p. We conclude that there is no advantage in increasing the density 
beyond the minimum in Fig. |5|a. 

3.3 a + b — i> c, homogeneous case 

An interesting situation arises when we turn off the reverse reaction. Then the problem 
consists of a- and 6-particles that diffuse and combine to form inert c-particles upon 
meeting. Equivalently, we may neglect species c altogether and think in terms of an 
annihilation reaction between 'particles' a and 'antiparticles' b. 

We assume the numbers of a- and 6-particles to be initially equal {pa{0) = pb{0) = p), 
which implies that they will remain so. Then, the density Pa{t) = Pb{t) obeys the rate 
equation dpa/dt = —kp'^, where k = Pi/t is the rate constant of the reaction a + b ^ c (of 
course P2 = 0). The solution of the rate equation is Pa{t) = p/{l + kpt), which behaves 
asymptotically in time as t~^. It is known, however, that, in microscopic simulations, the 
long-time behaviour of the system is determined by long-range density fiuctuations, which 
give rise to an asymptotic behaviour in two dimensions. Following Ref. [20|, this 



can be understood as follows: Initially = 0), particles a and b are distributed randomly 
with uniform probability. At a much later time t, particles have spread over a length 
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scale £o = (Dt)^/"^ due to diffusion and fluctuations on a longer length scale have not had 
time to dissolve. In a two-dimensional region of linear size Id, there are initially on the 
order of p£j^ particles of each species; however, due to statistical fluctuations, there will 
be deviations from this number. Assuming variations of one standard deviation, one can 
expect the difference between the numbers of a- and 6-particles to be of order (pi'jj)^^'^- 
Therefore, in this region of volume ~ i'^ , there is an excess density of a- or 6-particles on 
the order of p^^'^/io- For a given rate constant k, we can always choose the time t large 
enough so that the bulk of the particles will have annihilated and only the excess density 
will remain by that time. Thus, what remains of the initial system are regions occupied by 
one or the other species, with density of order £~^^ t^^/^; the reaction takes place only 
along the boundaries of these regions and is slower than under the well-mixed conditions 
assumed by the rate equation. Thus, the initial fluctuations in particle density give rise to 
a decay mechanism which is slower than t~^ and dominates at sufficiently long times. The 
rate equation, by neglecting fluctuations, does not account for this mechanism and predicts 
the wrong long-time behaviour. 

Simulations using our model reproduce asymptotically this result. At flnite times, we 
observe that the deviation from the result of the rate equations depends on the details of 
the experiment: (i) it decreases when the diffusion constant increases, since this limits the 
range of contributing fluctuations, and (ii) it increases with the rate constant, as the bulk 
of the particles annihilate faster and clusters develop sooner (Fig. ^, where r = 1). 

3.4 a + b ^ c, at an interface 

Moving away from homogeneous systems, we consider a system of a- and 6-particles 
which initially (t = 0) occupy the x > and x < halves of a two-dimensional {N^ x 
Ny) lattice respectively. For t > 0, they diffuse and react, producing c-particles (a + 
b ^ c); the latter diffuse and disintegrate into a- and 6-particles (c ^ a + 6). In the 
process, the a- and 6-particles diffuse into the regions of each other, while c-particles form 
a distribution peaked at re = 0. Here and in what follows, reaction rule II is used. We 
perform simulations, in which there are initially 50 000 particles on each half of a 500 x 500 
lattice (i.e. 0.4 particles/site) and other parameters are the same as before. Since the 
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Figure 6: Average density of a two-dimensional, homogeneous system reacting via a+b c 
according to rule II, as a function of time, for different values of the rate constant: simu- 
lation (solid curves) compared with rate equation (dotted curves). 



Figure 7: Width of reaction zone and of distribution of annihilation events for a 
two-dimensional, initially separated system reacting via a + b ^ c according to rule II, as a 
function of time: solution of reaction-transport equations (dotted, dashed and dot-dashed 
curves for 1, 4 and 10 d.p.r. respectively) compared with simulation (solid curves). 
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macroscopic problem is effectively one-dimensional for the particular initial condition, we 
define particle density as a function of time and the spatial variable x only, by averaging 
over the occupation number along the y-direction: 

= ^^A^„(x,y,t) . (3.2) 



'J/ 



y 



We wish to express quantitatively the size of the reaction front, i.e. the region around 
X = {) where reactions take place. This is obviously related to the spatial overlap of the 
densities of reacting particles. As a- and 5-particles spread on the lattice, their overlap 
steadily increases 0. Using pa{x,t), we define the width of the reaction zone, w, as the 
standard deviation of the product of the a- and 6-particle distributions: 

X{t) = —J2xpaix,t)pb{x,t) , w'^{t) = —Y,[x- X{t)f pa{x,t)pbix,t) . 

(3.3) 

The width obtained from these simulations is consistent, within statistical fiuctuations, 
with the asymptotic time dependence, namely t^/^, predicted by the reaction-diffusion 
equations. This behaviour is understood if we note that our choice of parameter values 
corresponds to a diffusion-limited regime. Then, particles essentially spread according to 
the t^/^ law expected from pure diffusion, with the relative abundance of different species 
determined locally by chemical equilibrium. If the reaction probabilities are much smaller 
than 1, then the time dependence of w displays more variety at times short on the time 



scale of the reactions, but the asymptotic behaviour remains the same [21 
3.5 a + b — > c, at an interface 

We now turn off the reverse reaction, as in the homogeneous case above, keeping all 
other parameters unchanged. In the diffusion-limited regime, a-b annihilation around 
X = proceeds so fast that new particles do not have the time to get there by diffusion. 
As a result, a depletion zone develops around x = 0, where the densities of a- and 6- 
particles are significantly smaller than their initial values. The width of this zone grows as 
can be shown from the reaction-diffusion equations that the width of the reaction 
zone, as defined previously, varies asymptotically with time as t^^^ |22|. For reasons that 



^"We always choose the size of the system so that the overlap does not reach the boundary during the 
simulation; the precise boundary condition is, therefore, irrelevant for the present discussion. 
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will become clear shortly, we define alternatively the width of the spatial distribution of 
annihilation events up to time t. Each time an a- and a 6-particle annihilate, we keep 
a record of the event and then we count at each lattice site all events that took place 
up to time t. The distribution of annihilation events is identical with the distribution of 
c-particles, provided the latter remain permanently at the site where they are produced. 
The new width, w' , is then defined by replacing in ( |3.3| ) the product of the a- and 6- 
densities by the density of the inert, immobile c-species. It can be easily shown that w' 
behaves asymptotically also as t^^^ according to the reaction-diffusion equations. 

In Fig. 1^ we compare the result of the simulation (solid curves) with that of the 
reaction-diffusion equations: the dotted curves describe the time evolution of the widths 
of the reaction zone {w, upper curve) and the distribution of reaction events {w', lower 
curve) for 1 d.p.r., while the dashed and dot-dashed curves describe the same quantities for 
4 and 10 d.p.r. respectively. The comparison is clearly facilitated by the better statistics 
in the case of w' . We notice that for 1 d.p.r. the result of the simulation lies above that of 
the reaction-diffusion equations; moreover, the former apparently grows with a higher time 
exponent than the expected 1/6. This is consistent with Ref. |1S]. It is clear that, with 



enhanced diffusion (4 and 10 d.p.r.), the result of the simulation converges systematically 
to that of the continuum equations. Continuing the simulation to greater times modifies 
slightly the time exponents, which depend further on the dimensionality of the lattice ||2^ 
and the reaction strength, but the following qualitative picture remains: in the absence of 
sufficient diffusive mixing, a and b annihilate less strongly and hence penetrate deeper into 
the regions of each other, resulting in a bigger additional widening of the reaction zone. 
This slowing down of the annihilation process is probably due to the inability of diffusion 
to destroy fluctuations beyond a certain length scale, as in the homogeneous case above. 
The analogy cannot be upheld, however, beyond the length scale set by the size of the 
depletion zone, which is absent in the homogeneous system. 

3.6 Complex reaction-diffusion systems 

We now consider reaction-diffusion systems which are significantly more complex than 
the ones considered so far. They involve species with different transport properties (e.g. 
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diffusion coefficients) and various chemical reactions, some of which are autocatalytic, i.e. 
they require the presence of a certain species in order to produce more of it. Autocatalytic 
reactions play an important role in biological processes. Systems subject to autocatalytic 
reactions can undergo phase transitions far from thermodynamic equilibrium. A system 
can be kept far from chemical equilibrium, for example, by suppressing reverse reactions 
and/or by an external supply of reactants. In such systems, the (unique) steady state, 
which is stable near equilibrium, may become unstable as certain parameters are varied. 
Then a phase transition to a new state may take place. Thus, an originally homogeneous 
state may become unstable and spatial concentration patterns (Turing structures) may 
develop spontaneously (see e.g. ^^). A crucial element in the formation of Turing 
structures is a significant difference in the diffusion coefficients of two species; the faster 
species, the inhibitor, hinders by chemical action the spreading of the slower one, the acti- 
vator, and the latter accumulates in restricted areas, creating a pattern of inhomogeneous 
concentration. The experimental observation of a sustained standing non-equilibrium 
chemical pattern has been reported recently The possibility of such distinct quali- 

tative behaviour makes autocatalytic systems attractive as a testing ground of the model 
developed here. It should be clear at the outset that the formation of structured states is 
predicted by the macroscopic diffusion-reaction equations. Linear stability analysis pro- 
vides, in fact, a critical value, dc, of the ratio of the inhibitor diffusion coefficient to that 
of the activator; when the ratio increases beyond dc, a range of Fourier components of the 
concentration become unstable and appear as spatial oscillations with the corresponding 
wavelengths. The details of the transition, such as the value of dc, may depend, however, 
on microscopic fluctuations and differences between the differential equation and CA ap- 
proaches may appear [27|. Our aim here is to show that our model is universal enough to 
describe qualitative phenomena of the kind described above. 



As a first example we consider the Brusselator |28], defined by the reaction scheme 



ki k2 fca 

A ^ X, B + X^Y + D, 2X + Y^2,X, X ^ E. (3.4) 

Here the concentrations of A and B are kept constant and D and E are continuously 
removed. In this case, X is the activator and Y the inhibitor. The Brusselator, albeit 
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relatively simple, displays the striking qualitative behaviour (concentration oscillations 
in time and space, nonlinear travelling waves) exemplified by the Belousov-Zhabotinski 
reaction |^^. We consider a one-dimensional system of length L. Following Ref. |25], 



one can put the reaction-diffusion equations for the unconstrained concentrations in a 
dimensionless form, 

^ = Dx^^+Ca-{Cb + 1)Cx + CICy (3.5a) 
^ = Dy'^ + CsCx-ClCy , (3.5b) 

by the transformation i = k^t, x = x/L, Dx,Y = Dxy /k^^L?' , Ca = {kiky^ /k^J'^)CA^ 
Cb = {k2/ki)CB and Cx,Y = (ks/k^y^'^Cxx ■ Since Cb is held fixed and D plays no active 
role, the second reaction in the original scheme can be effectively replaced hy X ^ Y, 
with rate constant /c2C'_B- 

Fig. ^ shows the result of a simulation of Eqs. (3.5) on a one-dimensional lattice of 
length I.OIL (1010 sites, A = X/L = 10^^). We choose the physicochemical parameters 



to have the values quoted in Fig. 7.13 of Ref. ||]: Ca = 2, Cb = 4.6, Dx = 1.6 x 10"^ 
and Dy = 3.75 Dx (dc = 3.05). The dimensionless concentrations Cx,y are computed by 
averaging the particle densities over cells of 10 sites and multiplying them by a scaling 
factor 7 ~ 0.2; to obtain Cx,y, we multiply by an additional factor (^4/^3)^/^. We 
constrain the boundary concentrations Cx and Cy to remain fixed at the values of the 
(unstable) homogeneous steady state, i.e. Ca and Cb/Ca respectively; this is done as 
follows: before every iteration step, we replace all particles in the first and last cells by 
uniform random distributions of particles having the necessary density. Initially X- and 
y-particles are distributed uniformly in the interior cell region of the lattice with densities 
of 9.2 and 10.7 particles/site respectively. The above value of 7 is chosen to give initial 
concentrations of 1.8 and 2.1 respectively. In the simulation we face a problem of competing 
reactions (e.g. a single X-particle may either convert to a y-particle or decay). Thus, in 
a naive implementation, in which all reactions are given a chance to take place one after 
the other, the order in which this is done would be significant. We solve the problem by 
dividing up reactions in groups of competing reactions and allowing randomly only one 
reaction from each group to take place during an iteration step (the probability of each 
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Figure 8: Simulation of one-dimensional Brusselator (solid curve) compared with solution 
of reaction-transport equations (dotted curve). 



Figure 9: Spatial concentration pattern obtained from simulation of two-dimensional 
Schnackenberg model; regions of high (low) density are indicated in red (blue). 
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reaction has to be multiplied of course by the number of members in its group). For the 
parameter values chosen, the initial homogeneous state is unstable and the system evolves 
to a structured state. The spatial structure shown is obtained after 800 000 iterations with 
time step f = k^r = 0.833 x 10~^. Obviously the statistics is still poor for a quantitative 
comparison with the continuum result, which, incidentally, has converged to a stationary 
value by the time shown here. 

The concentration pattern shown in Fig. |9| is produced with the reaction scheme 

fcl k2 ks 

A ^ X, 2X + Y ^ 3X, B ^ Y (3.6) 
fc_i 



(Schnackenberg model |3^). The simulation is performed on a 696 x 721 lattice of size 
1.10 X 1.14 (arbitrary units). Initially X- and y-particles are distributed randomly, with 
uniform density 1.55 and 0.59 particles/site respectively. Other parameters are chosen as 
in Fig. 2a of Ref. H: Cb = lAlkl'"^ /k^ky^ , Ca = 0.14 A;p//c_i/c2^^ Dx = lO-^/ciL^ 
(L being the linear size of the system). Finally we take Dy = 30 {dc — 20). The 
density of X-particles is shown here after 4 800 iterations, with r = 0.625 x 10^^ fc^"^. 
Time limitations have prevented us from running the simulation up to a time when the 
pattern becomes stationary. CA simulation of spatial patterns with the Schnackenberg 
model has recently been reported by another group p^. 



We have also performed simulations of the Se lkov reaction scheme, which was introduced in the context 
of glycolytic oscillations and is obtained from (|3.6|) by making all reactions reversible js^]. We obtain 
concentration patterns, in agreement with Ref. |33[ , both above and below the critical diffusion-coefficient 
ratio. Below the critical ratio, the homogeneous state is apparently destabilised by density fluctuations 
which effectively widen the range of unstable wavenumbers. 
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4 Conclusion and Outlook 



Reaction-transport processes were modelled in this work as a cellular automaton. Par- 
ticles are transported by executing a random walk on the sites of a regular lattice and are 
chemically transformed according to a local probabilistic rule. The microscopic random 
motion of the particles is manifested, at the macroscopic level, as a combination of advec- 
tion and diffusion. In particular, advection arises from a directional bias in the random 
walk, i.e. if particles have a relatively higher probability to move in the direction of the 
advection velocity. Chemical reactions are likewise modelled at the microscopic level: in 
the process of the random walk, we allow particles that meet at a lattice site the chance of 
'reacting', i.e. disappearing and leaving in their wake a set of new particles, the products 
of the reaction. The model is simple and general. The evolution equations are transcribed 
into simple computer code. This simplicity does not preclude, but, on the contrary, facil- 
itates the implementation of arbitrarily complex reactions and boundary conditions, in a 
physically transparent way. The results of the previous section demonstrate that the model 
presented here can successfully describe a wide variety of reaction-transport systems. 

In the continuum limit, the evolution equations of the discrete model go over to the 
standard reaction-transport PDE's, if certain conditions are fulfilled, namely if molecules 
of different species are uncorrelated (molecular chaos) and if particle density is smooth 
in space. These conditions arc, however, not imposed in our model; thus our simula- 
tions account for microscopic effects (e.g. fluctuations) that are typically averaged out by 
the continuum approach. The results of the discrete model were compared carefully to 
the solution of the reaction-transport PDE's in the case of a simple homogeneous system 
of particles. We encountered two sources of discrepancy: on the one hand correlations 
between the particles, that have no physical significance and are mere artifacts of the 
discretisation, and, on the other hand, statistical fluctuations that influence the long-time 
behaviour of reaction-diffusion systems. We emphasise the latter kind of discrepancy, 
which is of a fundamental nature and shows that microscopic fluctuations can influence 
qualitatively the evolution of macroscopic systems. The differences between the discrete 
simulation and the continuum approach in the time development of the reaction front 
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forming when the reacting species are initiahy separated is probably of similar origin. 
Macroscopic consequences of microscopic fluctuations were also seen in the case of auto- 
catalytic reaction-diffusion systems, near the threshold for the onset of non-equilibrium 
phase transitions leading to formation of spatially structured states. We have thus shown 
in concrete cases that the discrete model can he used to approximate systematically the 
respective PDE's, while, unlike the latter, it accounts for effects of explicitly microscopic 
origin. 

An important aspect of our approach is that the same macroscopic behaviour can be 
obtained with a variety of microscopic rules. This freedom is inherent in CA modelling, 
since the aim is to model only certain fundamental features of the microscopic world, but 
not the full detail of the dynamics. As examples, we gave two such rules: rule I, a simple 
rule leading to the standard rate law only in the limit of small densities, and rule II, a 
more complicated rule that, however, results in the standard rate law for any density. The 
computation time required by different rules in order to simulate the same macroscopic 
behaviour may vary. Our comparison of chemical rules I and II indicates that the latter 
affords better possibilities of optimisation, due to the wider range of possible particle 
densities. 

Finding the optimal microscopic rule is one way of enhancing computational efficiency, 
which is an essential task in practical applications. As a consequence of the unlimited 
number of particles per site, our model does not vectorise well when implemented on 
computers with vector architecture. In the present general form of the model, the best 
scalar performance is attained on a VAX-9000 and amounts to ~ 1 000 000 site updates 
per second (u.p.s.) for pure transport of one species, ~ 100000 u.p.s. for a + 6 ^ c and 
~ 20 000 u.p.s. for the Brusselator. This performance can be improved by varying the 
particle density, as we saw in Section ^. However, our model should be implemented on a 
massively parallel computer in order to achieve its full potential. The great advantage of 
our model is its ability to model a very wide range of reaction-transport problems. The 
performance figures quoted above are thus for a correspondingly general code. For specific 
applications, both the model and corresponding code can be appropriately tailored to 
optimise performance still further (see for example Ref. where 20 000 000 u.p.s. were 
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achieved on a CRAY-YMP for a specific reaction-diffusion process). 

An important feature of our CA approach is its inherent stabihty. Indeed our algorithm 
is stable even for values of r and A for which the corresponding FDE's, represented by 
Eqs. ( 2.34| ) and ( 2.33| ), are numerically unstable when solved by direct iteration with 
particle densities being floating point variables. We note that adding chemical reactions 
to a problem of pure transport can turn a stable FDE algorithm into an unstable one - 
for the same time and space discretisation 111. In such cases, where deterministic methods 
using floating point numbers fail, our stochastic approach, using integer variables, provides 
a stable method of solving Eqs. ( |2.34| ) and ( |2.33| ). The lack of general stability criteria for 
FDE in the presence of chemical reactions makes this guaranteed stability of our approach 
a valuable asset. 

Realistic applications will constitute the principal direction of further work on our 
model. These applications will be chosen on grounds of practical usefulness and so as to 
utilise the advantages of our approach: the guaranteed numerical stability and the capacity 
to treat arbitrary boundary conditions and chemical systems which are not necessarily in 
chemical equilibrium. Thus applications will include both problems at field scale, where 
differential equations are used widely to describe average quantities, as well at the scale 
of the pores, which form networks with very complex boundary conditions. As examples 
we mention the chemistry of mineral surfaces (precipitation/dissolution) and the effects 



of porosity changes on mass transport |34]. 

In conclusion, we have shown that the cellular automaton model we are proposing is 
capable of simulating simple as well as complex reaction-transport processes. We under- 
stand well its relation to other numerical methods for solving the corresponding differential 
equations and will proceed to apply it to realistic problems. Where other methods work 
as well, the relevant question will be that of the relative computational efficiency. More 
challenging will be, however, to identify areas of application where differential equation 
approaches are confronted with either serious numerical or even conceptual problems. 



^Stability is then achieved by refining further the time discretisation. 
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